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The coupled cubic nonlinear Schrodinger (CNLS) equations are used to study modulational insta- 
bilities of a pair of nonlinear ly interacting two-dimensional waves in deep water. It has been shown 
that the full dynamics of these interacting waves gives rise to localized large-amplitude wavepackets 
(wave focusing). In this short letter we attempt to verify this result numerically using a Fourier 
spectral method for the CNLS equations. 
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X> I. INTRODUCTION 

m 

Extremely large size waves (commonly known as freak, rogue or giant waves) are very common in the open sea 
or ocean and they pose major hazard to mariners. As early as 1976, Peregrine [10] suggested that in the region of 
oceans where there is a strong current present, freak waves can form when action is conncentrated by reflection into 
a caustic region. A variable current acts analogously to filamentation instability in laser-plasma interactions [7J [8] . 
Freak waves are very steep and is a nonlinear phenomena, hence they cannot be represented and described by a linear 
water wave theory. Zakharov [18] in 2009 has noted that in the last stage of their evolution, their steepness becomes 
'infinite', thereby forming a 'wall of water'. However, before such an instant in time, the steepness is higher than one 
for the limiting Stokes wave and before breaking the wave crest reaches three to four (sometimes even more) times 
higher than the crests of neighboring waves. The freak wave is preceded by a deep trough appearing as a 'hole in the 
sea'. On the other hand, a characteristic life time of a freak wave is short, typically ten of wave periods or so. For 
example, if the wave period is fifteen seconds, then their life time is just few minutes. Freak wave appears almost 
instantly from a relatively calm sea. It is, therefore, easy to appreciate that such peculiar features of freak waves 
cannot be explained by means of a linear theory. Even the focusing of ocean waves is a preconditions for formation 
of such waves. 

It is now quite common to associate appearance of freak waves with the modulation instability of Stokes waves. 
This instability (known as the Benjamin-Feir instability) was first discovered by Light hill [6] and the detail of theory 
was developed independently by Benjamin and Feir [1 and by Zakharov [14 . Zakharov showed slowly modulated 
weakly nonlinear Stokes wave can be described by nonlinear Shrodinger equation (NLSE) and that this equation is 
integrable [15] and is just the first term in the hierarchy of envelope equations describing packets of surface gravity 
waves. The second term in this hierarchy was calculated by Dysthe [3]. 

Since the pioneer work of Smith [12] , many researchers attempted (both theoretically or numerically) to explain the 
freak wave formation by NLSE. Among diverse results obtained by them there is one important common observation 
which has been made by all, and that is, nonlinear development of modulational instability leads to concentration of 
wave energy in a small spatial region. This marks the possibility for formation of freak wave. Modulation instability 
leads to decomposition of initially homogeneous Stokes wave into a system of envelope solitons, or more strictly 
quasi-solitons [T6J HZ] • This state can be called 'solitonic turbulence', or 'quasisolitonic turbulence'. 

In this letter, we consider the problem of a single solition in a homogeneous media, being subjected to modulational 
instability which eventually leads to formation of a system of soliton. We will show that the supercritical instability 
leads to maximum formation of soliton, concentrated in a small region. In going through subcritical instability 
the solitons coagulate to early stages of supercritical instability. Moreover, we investigate the full dynamics of 
nonlinearly interacting deep water waves subjected to modulational/filamentation instabilities, and we find that 
random perturbations can grow to form inherently nonlinear water wave structures, the so called freak waves, through 
the nonlinear interaction between two coupled water waves. The latter should be of interest for explaining recent 
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observations in water wave dynamics. 



II. NUMERICAL APPROACH 

In a pioneering work, a theory for the modulational instability of a pair of two-dimensional nonlinearly coupled 
water waves in deep water, as well as the formation and dynamics of localized freak wave packets was presented [9j [TT] . 
The two wave packets were investigated in the context of nonlinear optics by p~9j [20] , in Bise-Einsten condensates by 
[2T| [22] , in transmission by [23 , and in plasmas by many other authors [24U28] . 

Following [11], we consider the two-dimensional CNLS equations in the following form 

OA dA „ dA\ d 2 A D d 2 A d 2 A _ . 2| . ^.„ l2 , n 



and 



dt x dx y dy J a dx 2 dy 2 ^ dxdy 



C *^Z - C v^7. ) + + P^Ji ~ ^7^7. ~ t \ B2 \ B ~ 2 C \A\ 2 B = 0, (2) 



where A and B are the amplitudes of the slowly varying wave envelopes. The x and y components of the group 
velocity are given respectively by 

C x = ujk/2n 2 and C y = uj£/2k 2 
and the group velocity dispersion coefficients are 

a = u(2£ 2 - /c 2 )/8k; 4 , f3 = u(2k 2 -£ 2 )/8k 4 and 7 = -3^/4^ 4 . 
Also, the nonlinearity coefficients (as in [9 ) are given by £ = ujk 2 /2 and 

C = uj(k 5 - k 3 £ 2 - 3k£ 4 - 2k^K + 2fc 2 ^ 2 « + 2tn)j2K 2 (k - 2k). 
Here /c and £ are wavenumbers and uj is the wave frequency. They are related by uo = (the dispersion relation 



for deep water waves [5]) with g the acceleration due to gravity and n the wavenumber norm given by n = \Jk 2 + £ 2 . 
For detail description of the formulation of the problem, we refer to the original works [9j [TT] . 

The nonlinear strongly coupled system of equations (TJ and ([2| will be computed using a fast numerical algorithm 
based on the spectral method [2] [13] which is explained below. 



A. Fourier Spectral Method 

We first noticed that by letting S = A + B and D = A — B, the system ([I]) and Q becomes symmetric, obtained 
from and Q 

JdS „ dS „ dD\ d 2 S Q d 2 S d 2 D _ 
J 3D „ dD „ dS\ d 2 D B 2 D d 2 S ^ m 

where 

G(u, v) = l [(£ + 277) (|« + w| 2 + \u - v\ 2 ) u + (£ - 277) (|u + w| 2 - |u - w| 2 ) v] (5) 

O 

Then, we reduce the above system of PDEs Q-Q into a system of ODEs using the Fourier transform of u(x, y) 
which is defined by 

1 f CO [' 00 

F(u)(k x ,k y )=u(k x ,k y ) = — / e- i ^ x+k ^u(x,y)dxdy, (6) 

^ J — CO J — CO 
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with the corresponding inverse 

-j f'OQ f'OQ 

F- 1 (u)(x,y) = u(x,y) = — / e i{ - k * x+k ^u{k x ,k y ) dk x dk y . (7) 

The function u(k xi k y ) can be interpreted as the amplitude density of u for wavenumbers k xi k y . Now, we take the 
Fourier transform of both ([3| and Q as 



.dS t 
i—— 
dt 

.dD t 



(k x C x + ak 2 x + pk$) S-k y (C y + 7^) 5 = £>), 



(feC, + afc^ + /3A$ D — k y (C y + 7 fc x ) 5 = G(A 5), 
Letting fc^C^ + ak^ + /?fc^ = p and (C^ + 7^) = r ^ and ([9| can be written in the matrix form as 

d ( S \ ( p r 



dt\D 



r p 



S\ + G(S,D) 
D ) \G{d7s 



(8) 
(9) 

(10) 



Computing the eigenvalues and eigenvectors the solution to (10) can be written as 

S\ _ If e~ iXlt + e~ iX ^ - e ~ iXlt + e~ iX ^ \ ( 5(0) 
D J ~ 2 V -e~ iXlt + e" iA ^ e" iAlt + e" iA ^ / \ 5(0) 



(11) 



1 rt / e iX ^[G(S,D)-G(D,S)' 

2 7o I e iA2T \ G(S^D) + G(AS) 

with Ai = fc^C^ + afc^ + - fc^ (C^ + 7^) and A 2 = fc^C^ + afc^ + + fc^ (C^ + jk x ). 

We exploit the symmetry of the nonlinear function G from Q in developing a numerical procedure to solve the 



system of ODEs (10). 



B. Spatial discretization (Discrete Fourier Transform) 

We discretize the spatial domain Q = [— L/2, L/2] x [— L/2, L/2] into nxn uniformly spaced grid points Xij — {x\, yj) 
with Ax = Ay = L/n, n even, and L the length of the rectangular mesh Q. Given u(Xij) = Uij, i,j = 1, 2, • • • , n, we 
define the 2D Discrete Fourier transform (2DFT) of u as 



u kxky = AxAyJ2J2 e ~ iikxXi+kyVj)Ui ^ k x ,k y = -- + l, 
i=i 3=1 



and its inverse 2DFT as 



n/2 n/2 

^■ = (2^2 E E 1 

v 7 /c x = -n/2+l fc y = -n/2+l 



yi(k x Xi-\-kyyj ) 



i,j = 1,2, ••• ,n. 



(12) 



(13) 



In equation (12) and (13) the wavenumbers k x and k y , and the spatial indexes i and j, take only integer values. 



C. Temporal discretization 



We solve the initial value problem of the ODE system (10) using the classical fourth order Runge-Kutta (RK4) 
method and exact treatment for the linear part [2|. 

Given t max , we discretize the time domain [0, £ m ax] with equal time steps of size At with t n = 71 At, n = 0,l,2,---, 
and define S n = S(x,y;t n ) and D n = D(x,y;t n ). Initializing S n = S(t n ) and D n = D(t n ), we compute the Fourier 
transforms of the nonlinear terms T (g (t- 1 (S n ) , J 7 ' 1 (D 71 )^ and T [g fr- 1 n ) , J 7 ' 1 (S n ty , and advanced the 

ODE (10) in time with time step At using the explicit RK4 for the nonlinear part, together with an exact solution 
for the linear part as shown in (11). 
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D. Simulation setup 

The numerical code for the above procedure is implemented in FORTRAN 90 and executed on a Linux cluster. 
The initial profiles for A and B were taken as Gaussians, 



A(x,y;0) = (A + random(0 (lO' 3 / K))e- a{x2 + y2) 
B(x,y\0) = (B + random(O(10- 3 /K))e- aix2 + y2) 



(14) 
(15) 



In the simulations reported here, we used the parameter values #o = tt/6, g = 9.81, w = 0.56, k = 0.33, Ao = 0.1 /ft, 
Bq = Aq, 0, a = 1,0, L = 2 and a grid of 256 x 256 nodes in the computational domain [—1,1] x [—1,1] with the time 
step size At = 0.01. 

For each simulation we monitor the energies Qa{^) and Q#(£), calculated as 



/CO poo IL IL 

/ \A(x,y;t)\ 2 dxdy = Y,J2\ Ai i\ 2AxA y 
-co J-oo i=1 =1 



(16) 



/CO pOO IL 11 
/ \B(x,y;t)\ 2 dxdy = Y,Y,\ B ^ 2AxA y 
- CO J — CO -_i • i 



(17) 



i=l J=l 



Observing a finite energy will reveal stability of a solution. As soon as the solution becomes unstable, the energy 
diverges. When the solution dissipates the energy approaches to zero. 





FIG. 1: Energy evolution Qa{£) (left) and Qs(t) (right) in a typical 3000 min. simulation. 



III. RESULTS 



The results presented in this paper represent a preliminary study on dynamics of interacting nonlinear water waves. 
The problem considered here comprises the dynamics of nonlinear interacting water wave packets through solving the 
coupled system of equations ([T]) and ([2]). 

The results that are shown in Fig. 2 are all in dimensional units, where the two interacting waves initially have 
the amplitude A = B = 0.1/ft + ran, with ran representing a random low-amplitude noise, equal to 10 -3 /ft, in order 
to enhance instability. The results shown represent different time steps (starting on the left-hand panel and going 
downwards) for t = 300/cj, t = 600/cj, t = 900/cj then (right-hand panel) t = 1200/cj, t = 1500/cj; the last figure 
on the right-hand panel is at the same time as that above it but plotted from a different prospective reflecting the 
maximum growth rate in the y direction. For our simulations we have taken typical data from ocean waves [4 . The 
waves A and B in Fig. 2, then have the initial amplitudes \A\ = \B\ = 0.1/ n ~ 3 meters. From these figures, we 
see at t = 1500/cj (« 2680 seconds) that wave A focuses as a localized wave packets with a maximum amplitude of 
« 0.35/ft w 10 meters. We remark for considerable period after the initial step, waves A and B are qualitatively the 
same (with \A\ > \B\) before the nonlinear wave- wave interactions set in which results to wave break-up. 

In summary, we presented a numerical procedure to solve CNLS equations describing modulational instabilities of 
a pair of nonlinearly interacting two-dimensional waves in deep water. The simulation results of the full dynamical 
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FIG. 2: The interaction between two waves, with equal initial amplitudes \A\ = \B\ = 0.1 ft -1 which are propagating at an 
angle of — tt/6. A low-amplitude noise equal to 10 -3 / k is added to the initial amplitude in order to enhance the modulation 
instability. 

system reveals that even waves that are separately modulationally stable can, when nonlinear interactions are taken 
into account, give rise to novel behavior such as the formation of large- amplitude coherent wave packets with am- 
plitudes several times the initial waves. This behavior is quite different from that of a single wave (the case for the 
original Benjamin- Feir instability) which disintegrates into a wide spectrum of waves. These results are relevant to 
the nonlinear instability arising from colliding water waves thereby producing large-amplitude oceanic freak waves. 
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